#include <math.h>
#include <stdio.h>
#include <stdlib.h>

bool calc(const double s0, const double m1, const double m2, const double m3, const double m4, const double deg)
{
  /*
    1+2 -> 3+4
    2 : target. Scattering angle of 3 is theta (rad).
   */

  const double pi = 3.1415926536;
  const double eps = 1e-10;

  const double ss = s0*s0;
  const double theta = deg*pi/180;

  printf("\n");
  printf("sqrt(s): %f\n", s0);
  printf("m1: %f\n", m1);
  printf("m2: %f\n", m2);
  printf("m3: %f\n", m3);
  printf("m4: %f\n", m4);
  printf("theta: %f deg = %f rad\n", deg, theta);
  printf("\n");

  const double e1 = (ss-(m1*m1+m2*m2))/2/m2;
  const double p1 = sqrt(e1*e1-m1*m1);

  printf("e1: %20.6f p1: %20.6f\n", e1, p1);
  printf("\n");

  //================================
  const double p_cms = sqrt( (ss*ss + m1*m1*m1*m1 + m2*m2*m2*m2 -2*(ss*m1*m1+ss*m2*m2+m1*m1*m2*m2))/4/ss );
  printf("p_cms: %20.6f\n", p_cms);
  printf("\n");

  //=============

  const double AA = m1*m1 + m2*m2 + m3*m3 - m4*m4 + 2*m2*e1;
  const double BB = 2*m2 + 2*e1;
  const double CC = 2*p1*cos(theta);
  const double delta = AA*AA - (BB*BB-CC*CC)*m3*m3;
  const double denor = BB*BB-CC*CC;
  if(delta<0 || fabs(denor)<eps ){
    printf("delta<0 or denor<eps !! %f %f -- %f %f %f -- %f %f\n", e1, p1, AA, BB, CC, delta, denor);
    return false;
  }

  printf("\n----------------------------------\n\n");

  const double e3[2] = {(AA*BB + CC*sqrt(delta))/denor, (AA*BB - CC*sqrt(delta))/denor};
  for(int ii=0; ii<2; ii++){
    const double positivecheck = (BB*e3[ii]-AA)/CC;
    printf("(BE-A)/C = sqrt(E^2-m^2) %s 0!! A: %f, B: %f, C: %f, E: %f, BE-A: %f\n", positivecheck > 0 ? ">" : "<", AA, BB, CC, e3[ii], BB*e3[ii]-AA);
    if(positivecheck<0){
      printf("\n");     
      continue;
    }
    //----------
    const double p3 = sqrt(e3[ii]*e3[ii]-m3*m3);
    const double sintheta = p3*sin(theta)/p_cms;
    const bool kpass = (sintheta<=1);
    printf("%d: e3: %20.6f p3: %20.6f -- %s\n", ii, e3[ii], p3, kpass? "good" : "bad");
    if(!kpass){
      printf("\n");
      continue;
    }

    const double theta_cms = asin(sintheta);
    const double deg_cms = theta_cms/pi*180;
    printf("%d: theta_cms: %20.6f deg = %20.6f rad\n", ii, deg_cms, theta_cms);
    printf("\n");

    const double bg = p_cms/m2;
    const double beta = bg/sqrt(1+bg*bg);
    const double gamma = bg/beta;

    const double pl = gamma*p3*cos(theta)-bg*e3[ii];
    const double costheta = pl/p_cms;
    if( (costheta*costheta+sintheta*sintheta-1)>eps ){
      printf("thetacms not consistent!! %f %f -- %f\n", costheta, sintheta, costheta*costheta+sintheta*sintheta-1);
      return false;
    }
    else{
      printf("theta cms consistent from sin and cos!!!\n");
    }
    printf("\n----------------------------------\n\n");
  }
  
  return true;
}

int main(const int argc, char * argv[])
{
  double argf[argc];
  for(int ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
    argf[ii] = atof( argv[ii] );
  }

  if(argc!=7){
    printf("argc!=7!!\n");
    return 1;
  }

  calc(argf[1], argf[2], argf[3], argf[4], argf[5], argf[6]);

  return 0;
}
